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ABSTRACT 

We  consider  the  recent  algorithms  for  computing  fixed  points  or  zeros  of 

n> 

continuous  functions  from  R to  itself  that  are  based  on  tracing  piecewise- 
linear  paths  in  triangulations.  We  investigate  the  possible  savings  that  arise 
when  these  fixed-point  algorithms  with  their  usual  tri angulations  are  applied 
to  computing  zeros  of  functions  f with  special  structure:  f is  either  linear 
in  certain  variables,  separable,  or  has  Jacobian  with  small  bandwidth.  Each  of 
these  structures  leads  to  a property  we  call  modularity;  the  algorithmic  path 
within  a sinplex  can  be  continued  into  an  adjacent  simplex  without  a function 
evaluation  or  linear  programming  pivot.  Modularity  also  arises  without  any 
special  structure  on  f from  the  linearity  of  the  function  that  is  deformed  to 

n the  case  that  f is  separable  we  show  that  the  path  generated  by 

2 

Kojima’s  algorithm  with  the  homotopy  H coincides  with  the  path  generated  by 
the  standard  restart  algorithm  of  Merrill  when  the  usual  triangulations  are  em- 
ployed. The  extra  function  evaluations  and  linear  programming  steps  required 
by  the  standard  algorithm  can  be  avoided  by  exploiting  modularity. 
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SIGNIFICANCE  AND  EXPLANATION 


--V 


Problems  of  solving  systems  of  equations  in  several  variables  arise  fre- 
quently in  applied  mathematics  - in  solving  discretized  versions  of  boundary 
value  problems,  in  optimization  and  in  economics.  A class  of  algorithms  for 
such  problems  has  been  developed  in  recent  years  based  on  tracing  piecewise- 
linear  paths  in  a triangulation  of  the  domain  of  interest.  These  algorithms 
are  generally  called  fixed-point  algorithms  because  they  guarantee  convergence 
when  seeking  a fixed  point  of  a continuous  function  from  a compact  convex  set 
into  itself.  Besides  strong  global  convergence  properties,  these  algorithms 
possess  quadratic  convergence  to  a solution  given  sufficient  smoothness  condi- 
tions . 

The  main  drawback  of  these  algorithms  is  the  large  number  of  function 
. evaluations  required  compared  to  other  methods  that  may  not  be  as  robust.  Here 

we  alleviate  this  problem  by  showing  how  to  exploit  special  structure  in  the 
function  f whose  zero  is  sought  - either  linearity  in  certain  variables, 
separability , or  having  a Jacobian  matrix  whose  nonzero  entries  are  concentrated 
around  the  diagonal. 


lies  with  MRC,  and  not  with  the  author  of  this  report. 


EXPLOITING  STRUCTURE  IN  FIXED-POINT  COMPUTATION 


Michael  J.  Todd 


1.  Introduction. 

We  consider  fixed-point  algorithms  [2,3,4,5,10,15,18]  for  computing  a zero  of  a contin- 
uous function  f : Rn  -*•  Rn.  Such  algorithms  deform  a simple  function  f°  : Rn  ■*  Rn  into  f 
or  into  a piecewise-linear  approximation  to  f and  trace  the  zeroes  of  the  resulting  homotopy. 
Our  aim  is  to  see  how  such  algorithms  can  exploit  special  structure  in  f and  f°. 

While  our  approach  is  valid  for  other  algorithms,  in  particular  the  homotopy  method  of 
Eaves  [2]  and  Eaves  and  Saigal  [4]  and  the  algorithm  of  Garcia  [5] , we  shall  confine  ourselves 
to  the  restart  method  of  Merrill  [10].  Here  f°  is  an  affine  function  with  a unique  zero  at 
x°,  say.  Consider  the  homotopy  h : Rn  x [0,1]  -*■  Rn  defined  by  h(x,t)  - tf(x)  + (l-t)f°(x). 
Let  T be  a special  triangulation  of  Rn  x [0,1];  that  is,  all  vertices  of  T lie  in 
Rn  x {0,1}.  Next  let  i be  a piecewise-linear  approximation  to  h with  respect  to  T.  In 
other  words,  t agrees  with  h on  the  vertices  of  T and  is  affine  on  each  simplex  of  T. 

Merrill’s  algorithm  traces  the  piecewise-linear  path  of  zeroes  of  H.  starting  with  the 
known  zero  (x^,  0).  Under  certain  conditions,  this  path  ends  at  a point  (x\  1)  and  x^ 
is  thus  a zero  of  a piecewise-linear  approximation  to  f.  For  details  on  this  path  tracing, 
see  Merrill  [10]  or  Todd  [18] . One  can  then  restart  the  algorithm  with  a special  triangulation 
T'  of  smaller  mesh  and  a new  function  f°  whose  unique  zero  is  x1. 

Suppose  the  algorithmic  path  meets  an  (n+1) -simplex  a = <y°,...,yn+H  of  T.  Then 
within  this  simplex  the  path  is  linear.  Assume  that  as  the  path  is  traced  one  encounters  the 

— n 

face  of  a opposite  y . The  algorithm  requires  that  we  find  the  (n+1) -simplex  o'  of  T 
on  the  other  side  of  this  face,  i.e.  o’  = <y®,...,y^  y^ , , . . . ,yI1+*)  e T,  o'  * a.  We 

next  compute  My^)  and  hence  the  algorithmic  path'within  o'.  Our  interest  is  in  the  case 
where  t is  affine  in  o u o';  hence  l(y^)  and  the  algorithmic  path  in  o'  can  be  pre- 


1 


l 


dieted  from  information  known  at  a - indeed  the  path  in  o'  is  merely  an  extension  of  the 


straight  line  path  in  o.  In  this  ease  we  say  l is  modular  at  with  respect  to  a. 

In  section  2 we  describe  two  triangulations  of  Rn  x [0,1]  commonly  used  in  fixed-point 
computation,  for  which  modularity  is  easily  exploited.  In  terms  of  these  triangulations  the 
definition  of  modularity  takes  an  especially  simple  form  reminiscent  of  the  definition  of  modu- 
larity for  functions  on  a lattice.  In  section  3 we  show  how  modularity  follows  from  three 
structures  on  f:  linearity  in  certain  variables , separability,  or  the  property  of  having  a 
Jacobian  with  small  bandwidth.  Section  4 describes  how  we  exploit  modularity  in  the  restart 
algorithm. 

Sections  5,  6 and  7 discuss  further  the  special  structures  above.  In  section  5 we  show 

how  linearity  arises  naturally  in  the  zero  finding  problem  (or,  more  generally,  complementarity 

problem)  arising  from  a nonlinear  programming  problem.  Section  6 contains  the  surprising  re- 

2 

suit  that  Kojima's  algorithm  [8]  for  the  separable  case  using  his  homotopy  H generates  a 
path  identical  to  that  of  the  standard  restart  algorithm.  The  extra  linear  programing  pivots 
and  function  evaluations  required  by  the  latter  may  all  be  avoided  by  exploiting  modularity. 
Finally,  section  7 discusses  how  function  evaluations  can  be  reduced  if  the  Jacobian  of  f has 
small  bandwidth  - our  approach  here  follows  Curtis,  Powell  and  Reid  [1] . 


I 


! 
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2.  Triangulations. 

For  simplicity  we  confine  our  attention  to  triangulations  with  grid  size  1.  When  used  in 
fixed-point  computation,  these  triangulations  are  generally  scaled  and  translated  in  their 
first  n coordinates.  Identify  Rn  » [0,1]  with  {x  e Rn+1|o  <_  xn+1  i.  Removal  of  the  bar 
from  a vector  in  Rn  x [0,1J  denotes  its  projection  into  Rn.  Let  u1,...,un+1  be  the  unit 
vectors  in  Rn+1,  so  that  u* , . . . , un  are  the  unit  vectors  in  Rn. 

Now  let  y £ Rn  « {0}  have  y^  an  integer  for  1 <.  i <_  n.  Let  i be  a permutation  of 
{1, 2, . . . ,n+l) . Then  k^fy.n)  denotes  the  simplex  (y° , . . . .y"*1!  , where 


-0 

y 


- -i  — i-1 

y ; y - y 


r(i) 


1 < i < n+1 


(1) 


The  special  triangulation  is  the  set  of  all  such  k^fy^l’s  [18,  chapter  III). 

Next  let  y c Rn  x {0}  have  y an  odd  integer  for  1 < i < n.  Let  n be  a permutation 
of  {1,2,. . . ,n+l}  and  let  s £ Rn  x {1}  be  a sign  vector  - s^  = ±1  for  1 £ i <_  n.  Then 
j1<y#*»s)  denotes  the  simplex  <y°, . . . ,yn+1)  , where 


-0 

y 


— —i 

y > y 


— i-l  , - 
y + s 


i(i) 


tr(i) 


1 < i < n+1. 


(2) 


The  special  triangulation  is  the  set  of  all  such  j^fy.n.s) 's  [18,  chapter  Till. 

The  property  that  allows  us  to  efficiently  exploit  structure  in  these  triangulations  is 
that  each  vertex  of  a simplex  differs  from  its  predecessor  in  just  one  coordinate.  We  always 
assume  the  vertices  of  a simplex  in  or  are  ordered  as  above  to  simplify  the  use  of 
this  property. 


Let 

,-0  -n+1. 

o = <y  ,...,y  > 

E T,  T =■  Kj  or  J^,  and 

suppose  that  we  wish  to 

know  y3  so  that 

when 

yj 

— n 

replaces  yJ  in 

a,  another  sinplex  in  T 

results.  (We  say  y3 

is  the  replace- 

ment 

for 

— -j  - 

y in  o.)  Consider  first  K^.  We  have 

o 

■ 

•n 

*0  -n+1  -0  1 . 

.y  »y  — y +y,o' 

.-1  -n+1  -0. 

- <y  ,...,y  , y~ 

0 < j < n+1 

•j  -j-1  -)  . -j+1 

! y " y ~ y + y , a 

i 0 -j  — n+L 

i*  «*  <y  ....  ,y  , . . . ,y  ~ 

(3) 

j » n+1  i 

-n+1  — n -n+1  —0  , 

! y « y - y + y , ® 

, - n+1  —o  —tv 

■ <y  , y ,. . . ,y  > 
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In  each  case,  we  list  the  vertices  on  the  right  in  the  natural  ordering  given  by  (1) . 


Now  consider  J^.  We  have 


j = 0 s y 

0 < j < n+1  : y] 
*n+l 


-y°  + 2yl 


P'1  - p + ^>+1 


(4) 


n -n+1 
2y  - y 


j = n+1  s y 

0 -n=l. 

In  each  case,  the  natural  ordering  in  (2)  is  unchanged  - we  have  a - 'y  , . . . ,y  , . . . ,y  ' 

For  these  triangulations,  modularity  takes  an  especially  simple  form.  First  let  T = K^. 

— j — j— 1 — h — j+1 

Suppose  0 < j < n+1  and  let  h = v(j) , i = *(j+l) - Then  y = y + u , y 
p-1  + ^ and  y^  = p'1  + u1.  It  is  easy  to  see  that  l is  modular  at  p with  re- 

spect to  o if  n+1  ^ {h,i}  and  either  n 1(n+l)  < j and 


f(y^_1  + uh)  + f (y 


■3-1  + u1)  = f (y"*-1)  + f(y^_1  + uh  + u1) 


(5) 


or  it"1  (n+1)  > j and  (5)  holds  with  f°  replacing  f. 

The  similarity  of  (5)  to  the  condition  for  modularity  for  a function  defied  on  a lattice 
prompts  our  nomenclature.  (The  condition  given  above  is  almost  necessary  as  well  as  sufficient 
for  modularity.  Any  modularity  with  j € {0,n+l}  or  n+1  e (h,i)  can  only  be  due  to  a for- 
tuitous relationship  between  f and  f ^ . ) 

Now  suppose  T - j . If  0 < j < n+1,  a sufficient  condition  for  modularity  is  similar  to 
that  given  above  with  s^u^1  and  s^u1  replacing  u and  u in  (5)  . Now  suppose  j “ 0 
and  let  i = it ( 1)  . If  i = n+1  then  y°  is  the  only  vertex  of  a in  Rn  x {0}.  Thus  the 
end  point  (x1,!)  has  been  found  and  no  function  evaluation  or  linear  programming  pivot  is  re- 
quired.  Hence  suppose  i * n+1.  Then  i is  modular  at  y in  o iff 

f°(y°  - u1)  + f°(y°  + u1)  « 2f°(y°)  . (6) 

Finally,  if  j - n+1  let  h = 7r(n+l)  - since  y"+1  is  to  be  removed,  we  can  assume  h * n+1. 

Then  l is  modular  at  yn+1  in  a iff 


f(yn+1  - uh)  + f (yn+1+  uh)  - 2f(yn+1) 


(7) 


-4- 


Henceforth,  for  simplicity  in  notation  we  deal  only  with  the  case  that  T = K^.  The 
corresponding  results  for  J1  are  immediate.  As  we  shall  see,  there  are  advantages  to  3 - 

in  particular,  since  f°  is  affine,  (6)  always  holds.  From  now  on,  we  shall  also  abbreviate 
"t  is  modular  at  with  respect  to  o"  by  "y^  is  modular  in  o"  with  i understood. 


3.  Examples  of  Modularity. 


Here  we  show  how  modularity  follows  from  three  special  structures  on  f.  In  addition, 
modularity  arises  naturally  from  the  linearity  of  the  artificial  function  f°. 

3.1.  Linearity.  Identify  R^  » Rq  with  Rn,  where  p+q  * n,  and  suppose  f is  such 
that  f(v,w)  is  affine  in  w £ Rq  for  each  fixed  v e R^ . Section  5 gives  an  application  of 
such  a problem  to  nonlinear  programming.  We  then  have 

Lemma  3.1.  Suppose  o = (y°, . . . ,yI1+1>  = k^(y,t)  , and  f is  as  above.  Then  y^  is  modular 
in  a if  0 < j < it  ^(n+l)  - 1 or  if  tt  ^(n+l)  < j < n+1  and  both  h = ir(j)  and  i = n(j+l) 
are  greater  than  p. 

Proof.  If  j < tt  ^(n+l)  - 1 then  n+1  ^ {h,i}  and  it  ^(n+l)  > j.  Thus  we  only  need  to  check 

that  (5)  holds  with  f ^ replacing  f.  But  this  follows  from  the  linearity  (in  the  affine 

sense)  of  f°. 

If  j > it  *(n+l)  then  again  n+1  { {h,i}  and  we  must  verify  (5).  But  the  projections  of 
y^  y^  1 + uh,  y-*  1 + ui  and  y^  1 + uh  + u*  on  their  first  p coordinates  are  the  same. 
Hence  (5)  follows  from  the  linearity  of  f in  w for  fixed  v. 

Note  that  (6)  also  holds  if  f°  is  affine,  and  (7)  holds  if  h > p.  Hence  the  use  of 
allows  additional  exploitation  of  linearity. 

3.2.  Separability.  The  function  f is  called  separable  If  there  exist  functions 

f*  : R + Rn,  l<i<n,  such  that  f(x)  - ^ifi(x^) . Special  algorithms  to  compute  zeroes  of 

such  functions  have  been  studied  by  Kojima  18],  One  of  these  will  be  discussed  in  section  6. 

Lemma  3.2.  Suppose  a « <y° (....y*1*1)  = k^Cy.w)  and  f is  separable.  Then  y^  is 
modular  in  a if  0 < j < it  1(n+l)  - 1 or  it  1(n+l)  < j < n+1. 

Proof.  For  the  first  case  the  reasoning  follows  that  in  lemma  3.2.  Suppose 
* ^(n+l)  < j < n+1  and  let  h « n(j)  and  i ■ Tr(j+1)  . We  must  check  (5).  Consider  the 
equation 

f”<(yj-1  ♦ uh)  ) + <yj_1  + u1)  ) - fm(y^_1)  + f”((y1_1  + uh  + u1)  ) . 

In  m In  Iu 

For  m J (h,i)  all  four  terms  are  equal  and  the  equation  is  valid.  For  m *=  h,  the  first  and 


last  terms  are  equal  and  so  are  the  middle  two  - again  the  equation  is  valid.  Similarly,  it 
holds  for  m = i.  Summing  over  m gives  (5) . 

3.3.  Functions  with  small  bandwidth.  Suppose  the  pth  component  function  f^  of  f de- 
pends on  only  for  | p— q | < k.  Then  we  say  f has  bandwidth  2k- 1.  If  f is  continuously 

differentiable  and  its  Jacobian  matrix  has  bandwidth  m,  then  so  does  f. 

Lemma  3.3.  Suppose  a = <y^ , . . . ,yn+^>  = k^(y,ir)  and  f has  bandwidth  m.  Then  y^  is 
modular  in  a if  0 < j < it-1(n+l)  - 1 or  if  ir-1(n+l)  < j < n+1  and  | ir ( j)  - x(j+l)|  ^ m. 

Proof.  We  only  need  to  verify  (5)  for  | h— i | >_  m.  For  each  1 < p < n,  the  pth  com- 
ponent function  f is  affected  by  at  most  one  of  x,  and  x, . Hence  (5)  is  valid  for  f 

P li  i p 

since  the  terms  sure  equal  in  pairs.  Since  p was  arbitrary  the  lenina  is  proved. 

Section  7 discusses  further  how  the  occurrence  of  modularity  can  be  promoted  and  function 
evaluations  circumvented  when  f has  small  bandwidth. 


w 


4.  Exploiting  Modularity. 

Let  o « <y°,. . . ,yn+'H  be  a simplex  of  encountered  by  the  algorithm,  and  for  each 

j let  a^  e Rn+1  be  the  column  vector  (1,  tip)).  Then  a was  generated  because  some  face 

—a. 

of  a,  say  that  opposite  y , was  known  to  contain  a zero  of  £.  If  B is  the  matrix 

ak+1, .. . ,an+1l  , then  there  is  a solution  z > 0 to  Bz  = u1.  (In  fact,  to  cir- 
cumvent the  problems  of  degeneracy,  we  will  only  generate  this  face  if  B 1 has  lexicograph- 

ically nonnegative  rows.  For  simplicity,  however,  we  assume  nondegeneracy  throughout  - we 
suppose  B is  nonsingular  and  the  first  column  of  its  inverse  is  positive.) 

We  also  have  available  to  us  the  inverse  B-1.  The  analysis  below  assumes  that  B 1 is 
available  explicitly.  It  is  possible,  however,  to  assume  that  we  use  the  product  form  of  the 
inverse  and  still  exploit  modularity.1  A factored  form  of  the  inverse,  such  as  B = LU,  is 
not  suitable  for  the  use  of  our  ideas. 

— -1—1  — -lk 

The  standard  restart  algorithm  then  proceeds  as  follows.  Let  b = B u , a =■=  B a and 

find  j such  that  a.  > 0 and  b./a.  = min{b  /ala  > 0).  Here  and  throughout  this  section, 

3 j 3 r r1  r 

b and  a have  coordinates  indexed  0,...,k-l,  k+l,...,n+l,  to  correspond  to  the  indices  of 
the  columns  of  B.  Then  the  other  face  of  a containing  a zero  of  £ is  the  face  opposite 

— A 

y . Make  a linear  programming  pivot  step  to  obtain  B where  B is  the  basis  obtained  from 
B by  replacing  a^  with  ak.  Then  find  the  replacement  y"1  for  y^  in  o,  compute  £ (y  ^ ) 
and  introduce  a'’  -<1,  Up))  into  the  basis  B,  continuing  the  process. 

Now  suppose  that  ~p  is  modular  in  o.  Then  we  can  save  a function  evaluation  since 
. a^-1  - a^  + a^+1.  However  we  can  do  more  - we  need  not  make  the  pivot  step  to  obtain  B 1. 
Indeed,  suppose  first  that  j ^ {k-1,  k+1) . Let  B = la  , . . . ,a^ , . . . ,a  , a'  ,...,a  l be 
the  basis  obtained  from  B by  replacing  a1’  with  P . Then  the  inverse  of  B can  be  obtained 
from  B-1  as  follows.  The  (j+e)-th  row  of  B-1  is  the  (j+e)th  row  of  B 1 plus  the  jth 
row  of  B-1,  e - ±1.  The  jth  row  of  B_1  is  the  negative  of  the  jth  rcrw  of  B 1 . All 
these  relationships  follow  trivially  from  the  equation  a^  = a^  1 - a^  + a^+1.  Similar  changes 
are  made  to  a and  b.  The  total  work  involved  is  2n  + 6 additions  and  n+3  sign  changes. 


1I  am  grateful  to  Professor  Rcmesh  Saigal  for  this  suggestion. 
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Note  that  the  first  column  of  B X is  not  positive  - its  jth  entry  is  negative.  How- 

— — 1 ■ k 

ever,  since  a_.  was  positive,  it  is  now  negative;  if  we  express  u in  terms  of  B and  a , 

and  increase  the  weight  on  ak,  the  weight  on  a-'  increases  to  zero  and  beyond.  We  may  there- 

“■  If  ^ 

fore  simulate  the  entry  of  the  column  aJ  into  the  basis  B by  entering  a into  the  basis  B. 

We  therefore  recompute  the  minimum  ratio  min{br/a r/ar  > 0).  Notice  that  only  two  of 

these  quotients  (corresponding  to  r = j-1  and  r = j+1)  have  changed.  Suppose  the  minimum 
• ~ n * 

ratio  occurs  for  the  index  j ' . If  y is  not  modular  in  the  new  simplex 

o'  = <y  , . . . ,yJ , . . . ,y  > we  perform  a normal  pivot  step,  replacing  aJ  with  a . If,  how- 
ever, )p  is  modular  in  o'  and  j'  | {k-1,  k+1}  we  may  perform  the  same  trick  again.  Hence 
several  pivot  steps  may  be  saved  (or  rather,  replaced  by  trivial  pivot  steps  requiring  minimal 
arithmetic)  if  several  modular  vertices  are  encountered. 

We  must  now  deal  with  the  case  where  j e {k-1,  k+1}.'*'  Suppose  j = k-1.  Then  B-*^  = 

-1  k-2  -1  k-1  , -1  k k-2  k-1  , - i J 

Ba  -Ba  +Ba=e  -e  +a,  where  e denotes  a unit  vector  with  coordinates 

indexed  0,...,k-l,  k+l,...,n+l  and  ith  coordinate  equal  to  one.  We  now  want  to  express  uX 

, , 0 k-2  *k-l  k n+1  , — ,0  k-2 

in  terms  of  a ,...,a  , a , a ,..., a , or  equivalently,  b in  terms  of  e ,...,e  , 


- ^ k-2  k-1  - k+1 

a + e - e , a,  e , 


Suppose  the  appropriate  weights  are  X _,...,X  Then  we 

0 n+i 


aiuk-i  + V + xi  * bi  • 1 4 {k'2-  k_1}  '• 


ak-2(Xk-l  + V + Xk-1  + \-2  = bk-2  1 


3k-l(Xk-l  + V ‘ Xk-1  = bk-l  * 


We  know  a nonnegative  solution  with  Xfc_1  = 0(1^  = b^_1/a^_1)  and  we  want  to  increase  Xfc  . 
Notice  that  a^_^  is  positive,  since  j = k-1  was  chosen  in  the  minimum  ratio  test.  Hence 
X = X^_^  + ^ increases  with  X^_^ , from  the  last  equation.  Rewriting  (8)  we  obtain 

a.X  + X.  = b.  , i ^ {k-2,  k-1)  ; 

(^k-2  + ak-l)X  + Xk-2  = <bk-2  + bk-l*  ’ (9) 

<Vi  - 11  x + xk  = bk-i  • 

XThe  rest  of  the  section  is  very  technical.  The  reader  may  skip  to  section  5 without  loss  of 
continuity. 
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We  now  use  (9)  to  perform  a minimum  ratio  test  to  determine  which  vertex  next  leaves  the 

-k-2 

current  simplex  o'.  If  this  leaving  vertex  is  y we  find  the  band  of  non-unit  columns 
in  the  representation  of  the  current  basis  in  terms  of  B widening.  Let  us  therefore  progress 

.-^4.2  — a A g 

to  the  general  case,  where  y ,...,y  have  been  replaced  (in  that  order)  by  y ,...,y  , 

— k-1  — k-2  — r *k-l  -k-2  -r 

and  y , y , . . . ,y  have  been  replaced  (in  that  order)  by  y , y , . . . ,y  . Of  course, 

when  any  such  y^  is  replaced,  we  assume  that  it  was  modular  in  the  current  simplex,  or  the 

process  terminates. 

We  therefore  want  to  express  b in  terms  of  the  columns  of 


[e 


r-1  — r-1 

,e  , a+e 


k-1 


- k-2  k-1  - - k+2  k+1 

, a+e  - e , a , a+e  - e , . . . 


— , s+1  k+1  s+1 

. . ,a+e  - e , e , 


n+1, 
,e  ] 


(10) 


If  the  corresponding  weights  are  X„,...,X  , ,,  we  have 

u n+i 


a (Xr+***+X  ) + Xi  = b^ , i 4 {r-1. 


,s+l}  ; 


a ,(X  +***+X  ) + X , + X = b , 
r-1  r s r-1  r r-1 


a.  ( X +• * *+X  ) + X 
1 r s l+l 


bi(  r <_  i < k-1  ; 


(ak-i  - 1>(V*+X*-1)  + Vi(Xk+,,,+xs)  = bk-i  ; lll> 

WV'+V  + (ak+i  • “‘W—V  =bk+i  ’ 

a . (X  +•  • *+X  ) + X . , = b.  , k+1  < i < s ; 

1 r s l-l  1 

a (X  + • • • +X  ) + X + X , = b 
s+1  r s s s+1  s+1 

Suppose  yr  is  the  new  vertex  and  the  weight  X^  is  to  be  increased.  Adding  equations 

indexed  r through  k-1,  we  get  (ar+* • •+a^_1) (Xr+* • *+Xg)  - X^  = (br+  +bk_^) . As  we  shall 

see,  the  fact  that  y1  left  the  simplex  implies  that  (ar+*-  •+a]c_1>  is  positive.  Hence  we 

cam  increase  X = X +•  • -+X  instead  of  X . Adding  the  equations  indexed  k+1  through  s 
r s r 


we  obtain 


(Vl+,"tasl(Xr V ‘ Xs  = (W*+bs> 
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with  row  and  column  j removed.  It  follows 


Hence  B = Bj  Pj£.  where  P*”  is  Pkrs 

that  B 1 ■ (P~^g)  1 B^,  and  it  only  remains  to  invert  Pk^g.  But  this  is  easily  done.  Let 

Q.  be  the  inverse  of  P,  ; we  have  Q,  explicitly  as 

krs  krs  *krs  F 


Then  (P^^)  * ~ Q~^  and  thus  B 1 = . Obtaining  B 1 from  B^1  requires  only 

(s-r+4) (n+1)  additions/subtractions. 

Case  1(b)  . j e {r+1, . . . ,k-l} . First  do  a standard  linear  programing  pivot  step  to  ob- 
tain B^1,  where  B^  = A^  1.  Then  an  analysis  identical  to  that  above  shows  that 
B 1 = 1 Bf1*  where  1 is  Qkrs  with  row  j and  column  j-1  removed. 

Case  1 (c) . j t {k+1, . . . ,s-l} : analogous  to  case  1(b). 

_ , ...  . . , . . .0  r-2  r-1  k-1  r k-1  k-2  k-1 

Case  1(d) . j = r-1.  Let  A^  = [a  , . . . ,a  , a - a , a - a ,...,a  - a , 

a*  ak, . . . ,an+*]  . Then,  with  B^  - a'*1,  B^1  Can  be  obta*ned  from  B 1 by  adding  the 

(r-1)  st  through  (k-2)nd  rows  to  the  (k-l)st.  Next  let  B^  = A~k  1 - we  obtain  B"1  from 

B^  by  a standard  linear  programming  pivot  step;  notice  that,  from  (12),  the  pivot  element 

(ar_1+* • •♦ak_1)  is  positive.  Now  from  (10)  we  have  B = B^  Dkrg,  where  Dkrg  is  the  matrix 
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with  inverse 


r-1 


k-r 


•-1  1 1 


s-k-1 


-1 

1 


!••••!  1 


--1 


Hence  we  easily  obtain  B 

Case  1(e) ■ j - k.  Let  A 


si- 


la 


k-1 


k-1  k+1 


k-1  k+2 

a , a , 


n+l, 

, a ) . Then , 


~k  -1  -1 

with  " *2  ' B1  °*5ta^ne<^  from  B by  adding  the  (k+l)st  row  to  the  (k-l)st.  Next, 


let  B„ 


-k-1 


Aj'"  ~ - since  a^_^  + ak+1  “ 1 is  positive  from  (12)  , this  is  a standard  linear 


programning  pivot  step.  Now  B = B E.  , where  E.  is  the  matrix 

Krs  Krs 
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< 


an+1).  With  = A~^ , b”1  is  obtained  by  adding  the  (k+2)nd  through  sth  rows  of  B 1 

-k+1  -1 

to  its  (k+1) st  row.  Then  a standard  linear  programming  pivot  step  gives  us  (Aj  ) . From 

this  it  is  easy  to  obtain  B-1;  the  argument  is  analogous  to  case  1(d) . 


Case  1(g) . j = s+1.  Analogous  to  that  above;  we  start  now  by  adding  the  (k+2)nd 
through  (s+l)st  rows  of  B * to  its  (k+l)st  row. 

We  have  finally  completed  the  analysis  of  the  nonmodular  case.  Fortunately,  the  case  re- 
maining is  much  simpler. 
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~ n * i 

Case  2.  y or  yJ  is  modular  in  the  current  simplex  o' 


Case  2 (a)  . j 4 (r-1, ,s+l}.  In  the  original  B,  replace  a3  by  a3  = a 


j-1 


a-’  + a-*+1j 


the  update  of  B is  trivial.  Correspondingly  update  b and  a and  return  to  the  minimum 
ratio  test  of  (12)  - note  that  only  two  ratios  have  changed. 

Case  2(b)  . j = r-1.  Note  that  (12)  then  implies  that  (ar  i+""+ajc  j}  is  positive, 
thus  justifying  our  claim  above  (12)  that  (a^f'+a^  ^)  was  positive  at  the  previous  step. 

In  this  case,  r is  decreased  by  1 and  we  return  to  the  minimum  ratio  test  of  (12) . Only 
two  ratios  (for  r-2  and  r)  need  be  calculated. 

Case  2(c) . j = s.  Then  s is  decreased  by  1 and  we  return  to  (12)  - only  the  ratios 
for  s-1  and  s+1  are  updated. 

Case  2(d) . j = s+1.  We  increase  s by  1 and  return  to  (12)  . The  ratios  for  s and 
s+2  are  recomputed. 

Case  2 (e)  . j e {r+1, . . . ,s-l) . The  new  column,  assuming  j < k,  is  a-’  = a^  1 + a^+1-  a3. 

Expressed  in  terms  of  the  basis  B,  it  is  B(  a + e-'  2 - e^  ^ + e''  - ek  1)  . It  might  be  possi- 

ble to  handle  this  case  and  derive  the  appropriate  counterpart  to  (12)  to  make  the  minimum 
ratio  test.  On  the  other  hand,  the  complications  of  case  1 suggest  that  we  avoid  these  diffi- 
culties. Certainly  if  case  2(e)  occurred  several  times  in  succession,  the  resulting  basis  would 
be  very  hard  to  deal  with.  Hence  we  suggest  ignoring  the  modularity  and  returning  to  the  appro- 

* i 

pnate  subcase  of  case  1.  At  least  one  pivot  has  already  been  saved,  and  we  know  a . 

We  now  give  an  estimate  of  the  amount  of  work  that  is  saved  by  exploiting  modularity.  We 
suppose  that  the  difficulties  of  case  2(e)  have  been  surmounted.  Our  estimate  is  based  on  the 
measure  of  directional  density  of  a triangulation,  studied  in  [17,  191.  Given  x,  d e Rn,  let 

$(t)  = (x+td,  t)  e Rn  x [0,1].  Let  T be  a triangulation  of  Rn  * [0,1]  and  o^,  o 2 € T. 

We  write  a ^ + c>2  if  $(t-n)  £ d^,  ♦(t+n)  e o2  for  small  positive  p,  and  if  the  vertex  of 
o2  not  in  lies  in  Rn  x (i).  Let  ni  = i + |{t  e (0,1)  | dj  + °2  for  solne  °i'  °2  £ I • 

Then  n^n^)  measures  the  number  of  function  evaluations  of  f°(f)  and  the  associated  number 
of  linear  programming  pivot  steps.  Let  (T ,d)  be  the  average  of  n^  when  x is  uniformly 


distributed  in  a cube  of  side  2,  for  T ■ K and  J^.  We  then  have 
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Theorem  4.1  [19] 


(a)  N0(K1#d)  - Ii(d"  + (d^l)*)  ♦ lt<i  jld.-d^ 

(b)  Njl^d)  =■  1 + ^(d*  + (di-l)')  + lt<i  j|di-dj| 

(c)  N0(J1,d)  = lt  i(|d.|  + (d.-l)+  + td.+l)-)  + JA<j  J(  |di+dj | + |d.-d.|) 

(d)  N^.d)  = 1 + l±  y(|d.l  + (d.-l)'  + <d.+l)+)  + XA<j  J(]d.+d.|  + |d.-d.|)  . 

Here,  for  any  real  X,  X+  = max{0,X}  and  X = X+  - X.  Using  the  same  arguments  we  can  prove 
the  result  below.  Of  course,  since  i)>(t)  is  straight  we  may  call  all  vertices  modular.  How- 
ever, when  we  say  that  so  many  pivots  can  be  avoided  by  exploiting  modularity,  we  mean  only 
those  which  are  guaranteed  to  be  modular  by  the  linearity  of  f°  or  the  special  structure  of 
f. 


Theorem  4.2 


(i)  By  using  the  linearity  of  f°,  we  can  avoid  o£  the  P£vots  counted  in 

(a)  above,  and  ^ jldjJ  + Ii<j  j(|dj+dj|  + ldi-djl'  o£  the  pivots  counted  in  (c)  above. 


(ii) 

If 

f 

is  ! 

pivots 

counted 

1 in 

ed  in 

(d) 

above. 

(iii) 

If 

f 

is  ! 

(iv) 

If 

f 

has 

the 


p<i<j 


(iv)  If  f has  bandwidth  m,  we  can  avoid  in  (b)  and  (d)  that  part  of  the  pivots  counted  in 


the  sums  for  which  | i— j 1 ^ m. 
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« 


5.  An  example  of  linearity. 

Consider  the  nonlinear  programming  problem  of  minimizing  6(v)  subject  to  g(v)  = o, 
where  0:1^  -*•  R and  g:RP  -*  R^  are  continuously  differentiable.  The  problem  of  finding  a 
stationary  point  of  the  Lagrangian  of  this  program  is  that  of  finding  a zero  of  f,  where 
f(v,w)  - (V0(v)  + Vg(v)w,  g(v)).  Here  w e R°-  - clearly  for  any  fixed  v,  f is  linear  in  w. 

The  inequality  problem  (g(v)  <_  0)  has  traditionally  been  formulated  for  application  of 
fixed-point  algorithms  as  a zero-finding  problem  for  a point-to-set  mapping  defined  on  RP. 

The  algorithms  converge  slowly  in  this  case.  If  g is  affine  then  the  methods  of  Kojima  [7] 
and  Todd  [ 20]  provide  fast  algorithms. 

We  are  here  concerned  with  the  case  where  g is  not  affine.  For  simplicity  we  consider 
only  the  equality-constrained  problem  above,  but  our  analysis  is  easily  extended  to  inequality 
or  mixed  constraints.  Two  papers  have  addressed  this  problem:  Kojima  [9]  and  Todd  [20].  Both 
propose  methods  that  use  triangulations  of  RP  x Rq  x [0,1]  (explicitly  or  implicitly)  in  which 
the  mesh  size  for  the  w variables  is  large  or  infinite  and  does  not  shrink  as  the  iterations 
progress.  The  reason,  of  course,  is  that  a piecewise-linear  approximation  to  f with  respect 
to  such  a triangulation  can  be  made  arbitrarily  close  to  f. 

Fixed-point  algorithms,  besides  yielding  approximate  fixed  points  or  zeros,  also  provide 
approximations  to  the  Jacobian  of  f [16,  13] . The  exploitation  of  this  information  is  the 
key  to  obtaining  quadratic  or  faster  convergence  - see  Saigal  [12]  and  Saigal  and  Todd  [14] . 
Unfortunately,  the  algorithms  proposed  in  [9,  20]  do  not  give  good  approximations  to  the  partial 
derivatives  of  f with  respect  to  the  v-variables.  Indeed,  the  approximation  to  3f(v,w)/3v.. 
is  the  difference  between  the  function  values  at  the  two  vertices  of  the  final  simplex  that 
differ  by  u^ . If  w is  the  value  of  w for  these  vertices,  what  is  obtained  is  an  approxi- 
mation to  3f  (v,w)/3v., , and  w may  be  far  from  the  approximate  value  of  w. 

We  may  obtain  good  approximations  to  the  Jacobian  of  f by  using  a fine  grid  size  for  the 
w variables  as  well  as  the  v variables.  The  disadvantage  using  the  standard  restart  method 
is  the  large  number  of  pivots  to  move  from  one  region  of  w-space  to  another.  In  effect,  one 
has  a general  (p+q) -dimensional  problem.  It  is  generally  felt  that  the  gain  in  smoothness 
ccnpared  to  the  standard  p-dimensional  formulation  is  not  worth  the  corresponding  increase  in 


dimension. 
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However,  the  analysis  of  section  4 shows  that  a fine  grid  size  can  be  employed  for  the 
w variables  and  most  of  the  work  of  moving  in  the  w-space  can  be  eliminated  by  exploiting 
modularity.  The  advantage  of  smoothness  may  now  compensate  for  the  increase  in  dimension.  In 
any  case,  the  results  of  [12,  14]  show  that  asymptotically  quadratic  convergence  can  be  achieved 
under  reasonable  conditions;  because  of  linearity  in  w,  asymptotically  each  iterate  requires 
only  p+1  evaluations  of  V9  and  Vg. 


6.  Separable  functions. 

Recall  that  f:Rn  ■+  Rn  is  separable  if  there  exist  fX:R  •*  Rn,  1 <_  i <_  n,  such  that 

r»  Q 

f(x)  * ^ f (x^ . Clearly  any  linear  function  such  as  f is  separable.  Kojima  [8]  has 

suggested  two  algorithms  for  computing  a zero  of  such  an  f - here  we  restrict  our  attention 

2 

to  that  using  the  homotopy  H . 

0 •»  n 

Using  f,  f and  the  triangulation  of  R x [0,1]  we  construct  the  piecewise-linear 

function  l as  in  the  introduction.  On  the  other  hand,  for  each  i,  1 <_  i <_  n,  we  can  use 
f , f01  and  the  triangulation  of  R1  x [0,1]  to  construct  iNx^.t).  Define  £(x,t)  = 


y.  ii(x. ,t> . 

*•1  1 


Lemma  6.1.  f 


l. 


Proof.  Suppose  l|x,t)  = a e R . Then  there  are  *o"”'Xn+l  and  a simplex 


<y^ , ....y™*1)  £ such  that  (x,t)  = ? and 


n+l 


r+J, * + f3 

L]=0  j 


yn+1  X 

J=0  j 


(13) 


X.  ^ 0,  j = 0,...  ,n+l 


where  f-'  = f(y^)  if  » 1 and  f^  » f°(y^)  otherwise.  Suppose  a = k,  (y,x).  Let 

n+l  l 


{i  |l  _<  i _<  n,  n 1(i)  < it-1  (n+l)  } and  1 “ {1, . . . ,n}\I.  . For  i £ I, , define  z°  = y?. 


2 -o  . , ^0  1 

‘i  " yi  + X'  ti  = ti 

-1„.  . , . -1 


Xii  » J{X  |w  (i)  <,  j < * (n+l)} 

0 1-02-0. 


'i'  “i 


•'ll 


> 0 and 

4 

- 1. 

Also  define  X^Q  » X^ | 0 

(n+l)} 

and 

Xi2 

■ £{X^|ir  X(n+1)  < j < n+l} 

o 

t 

O -H 

and 

tX  • t^  - 1.  Also  define 

r1d)} 

and 

Xi2 

- E(X  < j < n+l). 

*1'  i 

-1,,, 


J{X  j | 0 <_  j < w 1(n+l) } , 


-j  -o 

Yi 

t * A 


Consider  the  ith  component  of  (x,t)  - y3 • Because  y^  « y?  for  j < n 1(i), 

Xi 

ij 


-1, 


y^  - yi  ♦ 1 for  j >.  * (i) , we  obtain  xA  « zA  + X^  z^  + Xi2  zi 


iO  fci  + Xil  *1  * Xi2  *i 


Now  let  f 


f1^)  if  y^+1  - 1 and  f 


Similarly 

ij 


f0i(73) 
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« 


a,  and  collecting  terms  yields 


rn+l  r»f| 

otherwise.  Then  (13)  gives  ).  „ X.  ).  , f 

‘■]=0  j ‘•1=1 

r n,,  iO  , . il  , i2.  , iO  Oi  0.  i2  _i,  2.  il  ,i,  1. 

L (xi0g  + ^iig  Xi2g  * = a'  where  g = f (Zj^)  . 9 = f (zi)  and  g = f (zj  if 

i=l 

t1  = 1,  giX  = f0i(zX)  otherwise.  Together  with  the  trivial  facts  that  V?  . X.,  = 1,  X..  > 0, 
i l i-k=0  lk  ik  — 

this  last  equation  says  precisely  that  H(x,t)  = a. 

The  lemna  gives  immediately 

2 

Theorem  6.2.  The  path  generated  by  Kojima's  algorithm  using  homotopy  H and  triangula- 
tion ( J^)  for  each  component  coincides  with  the  path  generated  by  Merrill's  algorithm  using 

triangulation  K^fJ^) . 

Of  course,  many  more  pivots  may  be  required  in  the  latter  algorithm.  Indeed,  Kojima's 
algorithm  works  with  the  linear  system 

Zi-1  ll 0 \kgik  - 0 

to  \k‘i  - * 

^k=0  Xik  = 1 


(14) 


Xiki°  • 


During  a single  pivot  in  (14)  , it  is  possible  that  several  X^'s  for  i £ 1^  = 

{ i ' | X < 1 - t)  change  their  relative  magnitudes.  From  the  correspondence  between  X/s 
and  • s in  lenrna  6.1,  it  can  be  seen  that  each  such  change  requires  a change  in  the  permu- 
tation ir  and  hence  a pivot  in  (13)  . 

However,  the  path  of  either  algorithm  is  straight  in  a piece  corresponding  to  a single 
pivot  in  (14) ; hence  all  the  additional  pivots  correspond  to  modular  vertices  and  can  be  per- 
formed by  the  techniques  of  section  4.  With  the  exception  of  pivots  of  type  case  2(e)  of 
section  4,  we  can  therefore  simulate  Kojima's  algorithm  in  a standard  restart  algorithm  - note 

that  the  latter  requires  only  a linear  system  of  order  n+1  rather  than  3n  as  in  (14) . Of 

course,  n of  the  constraints  of  (14)  are  of  generalized  upper  bound  form. 

The  theorem  also  inplies  that  the  asymptotic  results  of  quadratic  convergence  in  Saigal 
[12],  Saigal  and  Todd  [14]  apply  to  Kojima's  algorithm  also.  No  such  results  are  available 
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for  the  homotopy  H3  18]  - the  possibility  of  acceleration  seems  a powerful  reason  for  choosing 
H2  over  H3. 

A slight  modification  is  suggested  for  accelerating  in  the  separable  case.  In  the  standard 

algorithm,  having  obtained  an  approximate  zero  x3 , we  choose  f*3  to  have  a zero  at  x3 . We 

also  translate  the  triangulation  so  that  (x3,0)  is  in  the  center  of  the  face  of  the  starting 

simplex  <v° , ....y"*3)  opposite  y*1*3.  This  is  done  by  arranging  that  x3  = -r-(y°  + y11)  + — 

2n  n 

y"_3  y3.  Then  if  all  zeroes  of  l have  projections  within  e/2n  (in  the  norm)  of  x3, 

2 

only  n+1  function  evaluations  and  linear  programming  pivots  are  required  to  obtain  x ; here 
e is  the  grid  size  of  the  triangulation  or  used.  With  f separable,  it  is  prefera- 

ble to  translate  the  triangulation  so  that  x3  = (i  - 6) (y°  + yn)  + [26/(n-l) ] £"”3  y3 , where 
6 is  small  and  positive.  In  this  case,  as  long  as  all  zeroes  of  i have  projections  within 
(4  - 5) e of  x3,  only  n+1  function  evaluations  and  "nonmodular"  linear  programming  pivot 
steps  are  required.  Many  more  pivot  steps  may  be  needed,  but  they  will  all  be  of  the  trivial 
type  considered  in  section  4.  Thus  the  algorithms  can  be  accelerated  more  safely  in  the 
separable  case. 
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7.  Functions  with  small  bandwidth. 


Recall  that  f has  bandwidth  m » 2k-l  if  each  component  function  f^  depends  on 
only  for  | i— j | < k.  We  assume  6k-4  n.  Functions  with  small  bandwidth  arise  naturally  in 

discretizations  of  boundary  value  problems  (see  chapter  1 of  [11]).  In  addition,  if  f is 
"sparse"  (i.e.  its  Jacobian  matrix  is)  we  may  permute  the  coordinates  of  the  domain  and  range 
space  to  attempt  to  obtain  a small  bandwidth.  See  [6]  for  computational  complexity  results 
for  minimizing  the  bandwidth  by  such  permutations. 

Our  approach  is  based  on  that  of  Curtis,  Reid  and  Powell  [1].  Our  aim  is  not  just  to  avoid 
function  evaluations  but  to  promote  the  occurence  of  modularity.  Fortunately  these  goals 
suggest  the  same  strategy. 

We  would  like  to  encounter  simplices  k^ty.x)  such  that  adjacent  members  of  x differ  by 
at  least  m.  While  the  algorithm  is  deterministic,  we  can  at  least  control  the  starting  simplex. 
We  choose  this  to  have  associated  permutation  x = (1,  m+1,  2m+l,...,2,  m+2,  2m+2, ....... ,m, 

2m,...,n+l).  Then  the  first  step  of  the  algorithm  that  does  not  transfer  a vertex  from 
Rn  x (0)  to  Rn  x {1}  or  vice  versa  entails  a modular  pivot.  Indeed,  since  each  simplex 
change  affects  the  permutation  only  by  an  adjacent  transposition,  it  is  likely  that  several 
modular  pivots  will  be  encountered. 

Under  reasonable  smoothness  conditions  on  f,  we  know  that  asymptotically  only  n+1  sim- 
plices will  be  encountered.  These  all  correspond  to  k^(y,x),  with  x containing  the  indices 
1 through  n in  the  same  order  as  in  x,  and  n+1  moving  from  the  last  position  to  the  first. 
In  this  case  one  can  also  economize  on  function  evaluations. 

Of  course,  if  scalar  function  evaluations  are  much  cheaper  than  vector  evaluations  one  may 

—4  ,4.4 

economize  at  each  step.  If  y is  to  leave  the  simplex  and  be  replaced  by  y , then  f(y  ) 

can  be  obtained  from  f(y^  *)  or  f(y^+S  by  making  only  the  m scalar  function  evaluations 

in  which  it  is  know  to  differ.  Frequently,  however,  economies  in  common  expressions  and  the 

need  to  control  subroutine  calls  suggest  that  vector  function  evaluations  are  more  efficient. 

In  this  case  one  may  proceed  as  follows. 


Starting  with  <y^, . . . ,yT'+*'>  = k^(y,n)  we  first  calculate  f(yn+^).  If,  as  hoped,  y*1 
leaves  the  simplex,  it  is  replaced  by  yn;  we  need  f(yn+^  - u”*n*)  = f(yn  . Instead  we 
calculate  fly"  1)  . This  evaluation  gives  us  automatically  fty*  ) , f(y7'  *rl'+1),..., 

f(yn  ).  If,  as  we  again  hope,  yn  displaces  y”  1,  then  we  need  f(yn  2) ; but  this  is  al- 
ready known.  It  follows  that  if  the  sequence  of  simplices  is  as  hoped  (and  asymptotically 

. ,,  . , , . _ „ 0 ir_1(l)-l  t_1(2)-l  -1  , n+1 

guaranteed)  we  need  only  evaluate  f at  y = y , y , . . . ,y  and  y , 

a total  of  m+1  function  evaluations.  In  any  case,  we  have  not  made  any  extra  evaluations. 

It  is  clear  that,  at  a general  stage  of  the  algorithm,  we  may  also  try  to  guess  the  next 
few  vertices  that  will  be  generated  and  make  a function  evaluation  that  will  give  us  the  value 
of  the  function  at  each  of  these  vertices.  Also,  the  more  general  grouping  idea  of  Curtis, 
Powell  and  Reid  (1)  can  be  exploited  in  an  analogous  way  - asymptotically  only  m +1  function 
evaluations  are  necessary  each  cycle,  where  here  m is  the  number  of  groups. 
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